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In the last decade, minimal kinetic models, and primarily the Lattice Boltmann equation, have 
met with significant success in the simulation of complex hydrodynamic phenomena, ranging from 
slow flows in grossly irregular geometries to fully developed turbulence, to flows with dynamic 
phase transitions. Besides their practical value as efficient computational tools for the dynamics 
complex systems, these minimal models may also represent a new conceptual paradigm in modern 
computational statistical mechanics: instead of proceeding bottom-up from the underlying micro- 
dynamic systems, these minimal kinetic models are built top-down starting from the macroscopic 
target equations. This procedure can provide dramatic advantages, provided the essential physics 
is not lost along the way. For dissipative systems, one such essential requirement is the compliance 
with the Second Law of thermodynamics. 

In this Colloquium, we present a chronological survey of the main ideas behind the Lattice Boltz- 
mann method, with special focus on the role played by the H theorem in enforcing compliance 
of the method with macroscopic evolutionary constraints (Second Law) as well as in serving as a 
numerically stable computational tool for fluid flows and other dissipative systems out of equilib- 
rium. 
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The quest for a better understanding of the macro- 
scopic world in terms of underlying "fundamental" mi- 
croscopic laws is a timeless issue in the history of science 
and natural philosophy. Modern science has provided 
an admirably powerful theory, and mathematical tool as 
well, to address this issue in a sensible and productive 
way, a gift named Newtonian mechanics. Newtonian me- 
chanics is a theory of amazing depth and breadth, go- 
ing as it does from scales of planetary motion, all the 
way down to molecular trajectories, encompassing almost 
twenty orders of magnitude in the process. 

Besides an array of practical results, the application 
of Newtonian mechanics at the molecular scale generates 
a profound puzzle: the origin of irreversibility and the 
nature of time itself. Newton's equations are manifestly 
reversible, that is, invariant under time and velocity in- 
version, which means that molecular motion is basically 
like a movie which can be indifferently rolled forwards 
or backwards in time with no loss of information. This 
is in blatant contrast to our daily experience of an in- 
exorably one-sided nature of time evolution: the arrow 
of time travels only one-way. This puzzle has remained 
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outstanding for more than a century, and is still open 
to a large extent, and in any case is beyond the scope 
of the present work. Besides its profound philosophical 
implications, irreversibility also bears upon the practical 
question of predicting the macroscopic behaviour of com- 
plex systems comprising a huge number of (non-linearly) 
interacting individual units. The application of Newto- 
nian mechanics to the molecular scale inevitably leads 
to a gigantic wave of computational complexity due to 
the enormous number of atoms/molecules constituting 
macroscopic systems. This problem is circumvented by 
formulating continuum models, typically based on partial 
differential equations describing the space-time evolution 
of a few macroscopic fields, such as fluid density, pressure, 
temperature, and so on. This approach is indeed quite 
successful, but when confronted with complex systems 
out of equilibrium, e.g., fully turbulent flows, it shows 
clear limitations. The point is that, owing to strong non- 
linearities and multi-dimensionality, the aforementioned 
partial differential equations are often just too compli- 
cated to be solved by even the most powerful numeri- 
cal techniques. It therefore makes sense to go back to 
Newtonian-like dynamics, namely large sets of ordinary 
differential equations, and develop minimal fictitious par- 
ticle dynamics designed so as to relinquish as many mi- 
croscopic details as possible without corrupting the ulti- 
mate macroscopic target. 

In this work, we shall be concerned precisely with this 
type of modelling strategy. In particular, we shall turn 
our attention to alternative ways to gain understanding 
about the predictability of macroscopic phenomena out 
of hyperstylized "Newton-like" microscopic models. We 
hasten to add that these alternative routes are highly 
influenced by advances in computational modelling, and 
therefore they naturally fit into the general framework 
of computational statistical mechanics. Like their real- 
life physical counterparts, these hyperstylized models are 
required to display irreversible behaviour as a basic req- 
uisite of stability, whence the importance of designing 
them in compliance with the second law of thermody- 
namics. In this Colloquium we shall be concerned with 
the Lattice Boltzmann equation, a minimal form of the 
Boltzmann equation which retains just the least amount 
of kinetic information neeeded to recover correct hydro- 
dynamics as a macroscopic limit. The Lattice Boltzmann 
equation has proved quite effective in describing a vari- 
ety of complex flow situations using a very simple and 
elegant formalism built upon the aforementioned hyper- 
stylized approach. 



II. STATISTICAL MECHANICS BACKGROUND 

The theory of the Lattice Boltzmann equation belongs 
to the general framework of non-equilibrium statistical 
mechanics. In this section we shall therefore present a 
brief review of the main cornerstones of classical sta- 
tistical mechanics, starting from the most fundamental 



(atomistic) level, all the way up to the macroscopic level. 

For most practical purposes, our ability to predict the 
behaviour of the world around us depends upon the time 
evolution of macroscopic variables, e.g., pressure, tem- 
perature, flow speedsr, etc., which result from the col- 
lective average over of an enormous amount of individ- 
ual trajectories. Since we can only experience average 
quantities, it makes sense to think of mathematical for- 
mulations dealing directly with these average quantities: 
which is the chief task of Statistical Mechanics. 



A. The BBGKY hierarchy 

The traditional route to this end is the celebrated 
BBGKY (Bogoliubov-Born-Green-Kirkwood-Yvon) hier- 
archy (Cercignani, 1975, Liboff, 1998), leading from 
atomistic equations to fluid dynamic equations for the 
macroscopic variables, typically the Navier-Stokes equa- 
tion of fluid flow (Landau and Lifshitz, 1953). 

The BBGKY path is based on the four basic levels (see 
left branch of Figure 1): 

• Atomistic level (Newton-Hamilton) 

• Many-body kinetic level (Liouville) 

• One-body kinetic level (Boltzmann) 

• Macroscopic level (Navier-Stokes) 

Let us discuss these four levels is some more detail. 



B. The atomistic level 

The atomistic description of (classical) macroscopic 
systems is based on Newtonian mechanics. The math- 
ematical problem generated by Newtonian mechanics is 
to solve a set of N non-linear ordinary differential equa- 
tions: 



with the initial conditions: 

Xi(t = 0)=x Qi (2) 
Vi(t = 0) = v 0i , i = l,N 

where N is of the order of Avogadro's number Na ~ 
6 x 10 23 . In the above, mi are the molecular masses, Vi — 
dxi/dt the molecular speeds and f\ is the force acting 
upon the i-th molecule due to intermolecular interactions 
(Huang, 1987). 

The application of Newtonian mechanics to the molec- 
ular world prompts a daunting wave of computational 
complexity. A centimeter cube of an ordinary substance, 
say water, contains of the order of Avogadro's number of 
molecules. Keeping track of the motion of these many 
molecules in the way portrayed by Laplace, namely by 
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tracing in time the 6N variables Xi(t) and Vi(t) turns 
out to be a scientific chimera. Even assuming one has 
enough capacity to store so much information, one would 
still be left with the problem of dynamic instabilities in 
phase space: any tiny uncertainty in the initial micro- 
scopic state would blow up exponentially in time, thereby 
shrinking the predictability horizon of the system virtu- 
ally to zero. It is a great gift that such a huge amount of 
information, besides being unmanageable, is also needless 
as well, as we shall see in the next section. 

C. Many-body kinetic level 

The atomistic level deals with molecular positions 
and speeds and is governed by the Newton-Hamilton 
equations which describe a world of trajectories. The 
iV-body kinetic level deals with distribution functions 
/jv(xi, Vi . . .xn,vn, t), namely smooth fields describing 
the joint probability to find molecule 1 at position X\ 
with speed vi, and molecule 2 at position x 2 with speed 
V2, and so on up to molecule N around position x^ with 
speed vn, all at the same time t. Trajectories are here 
replaced by the notion of phase-space fluids obeying a 6iV 
dimensional continuity equation, known as the Liouville 
equation: 

N 

[d t + J2v i -d Xi +a i -d Vi ]f N = (3) 
i=i 

where = Fj/m, are the molecular accelerations. The 
underlying assumption is ergodicity: the time spent by 
the trajectory of the 6N dimensional coordinate P{t) = 
[a;i(i) . . . vtj{t)] in a given differential volume element Ar 
of phase space, is proportional to the measure of Ar. 

The Liouville equation does not by any means re- 
duce the amount of information to be handled via the 
Newtonian approach. In fact, since /jv is a continuum 
6N dimensional field, the amount of computational in- 
formation blows up exponentially! The Liouville equa- 
tion represents nonetheless a very valuable step, not 
because we can solve it, but because it sets the stage 
for a very elegant and powerful procedure to consis- 
tently eliminate irrelevant information. Simply inte- 
grate Jm over unwanted single-particle coordinates, to 
define low-order reduced distribution functions /m = 
fi2...M<N = J fi2...Ndz M +i—dz N , where dz k = dx k dv kl 
k = M + 1, . . . N. The result is a chain of equations: 

M 

[dt + ^v t ■ d Xt + ai ■ d Vi ]f M = C M (4) 

known as the BBGKY hierarchy. Note that the right- 
hand side collects the effects of intermolecular interac- 
tions. In the presence of a 6-body potential, Cm only 
involves b upper-lying distributions /m+i, •■•/m+6- For- 
tunately, most interesting macroscopic observables, such 
as density, pressure, temperature, and energy, only de- 
pend on 1 or 2-body distributions, so that our efforts can 



be channeled into the bottom floors, M = 1,2, of the 
BBGKY hierarchy. 

D. The Boltzmann equation and the ff-theorem 

The most important one-body equation is the cele- 
brated Boltzmann equation: 

d t f + v-d x f + a-d v f = C[f,f] (5) 

where f(x,v,t) is the probability density of finding a 
classical pointlike particle at position x at time t with 
speed v. The left hand side represents free streaming in 
phase space (x,v) and the right hand side denotes the 
effects of binary collisions, typically a very complicated 
integral operator encoding the details of molecular inter- 
actions. The Boltzmann equation relies on the famous 
molecular chaos ( "Stosszahlansatz") assumption: 

hi{xi,vuX2,V2,t) = f{x- L ,v- L ,t)f(x 2 ,v 2 ,t) (6) 

which asserts the absence of correlations between 
molecules entering a binary collision. It is precisely this 
arbitrary-if plausible- assumption which breaks time- 
reversal symmetry, since it is clear that after a col- 
lision molecules must be correlated because of mass- 
momentum- energy conservation. The essence of the 
molecular chaos assumption is that these post-collisional 
correlations decay exponentially fast in time so that the 
probability of the two particles to collide with each other 
again in a correlated state after any finite time lapse is 
virtually zero. Breaking time-reversal symmetry opens 
the door to irreversible behavior, and one of the most 
profound of Boltzmann's contributions to statistical me- 
chanics rests with his discovery of a quantitative measure 
of irreversibility: the celebrated _ff-theorem. This quanti- 
tative measure of irreversibility is provided by the Boltz- 
mann i7-function (in what follows, we term it also as the 
entropy function; the physical entropy being S = —k^,H): 

H{t)= J f(x,v,t)lnf(x,v,t)dvdx (7) 

which was shown by Boltzmann to be monotonically non- 
increasing in time, dH/dt < 0, regardless of the de- 
tails of the collision operator. The i7-theorem stands 
out as a conceptual bridge between micro- and macro- 
dynamics. Yet, it is difficult to think of a more debated 
and controversial issue in theoretical physics (see Wchrl, 
1978). We will not delve here into the details of the vari- 
ous paradoxes which were raised against the i?-theorem, 
nor shall we discuss the fact that Boltzmann derived his 
theorem without demonstrating under which conditions 
his equation, a complicated integro-differential initial- 
value problem, does indeed have solutions. While leav- 
ing mathematical rigor somehow behind, the _ff-theorem 
is nonetheless a monumental contribution to modern sci- 
ence, since it showed for the first time the way to a grand- 
unification of two fundamental and hitherto disconnected 
domains of science: Mechanics and Thermodynamics. 
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The practical importance of Boltzmann-like equations 
was only to be furthered by modern developments in 
theoretical physics, primarily the emergence of the fun- 
damental notion of quasi-particles as collective excita- 
tions of non-linear field theories (Kadanoff and Baym, 
1962). By shifting the focus from actual particles (real 
atoms or molecules) to quasi-particles, the Boltzmann 
equation goes way beyond the original framework from 
which it was derived (i.e., rarefied gas dynamics) and 
spreads nowadays its wings across a huge variety of 
fields in statistical mechanics, including neutron and ra- 
diation transport, electron transport in semiconductors, 
hadronic plasmas, and many others. Quasi-particles play 
a central role also in the top-down approach to statistical 
mechanics that is to be described shortly. 

E. The macroscopic level 

Macroscopic observables such as fluid mass density, 
speed and energy density are obtained from the one-body 
kinetic distribution by integration over velocity space: 

p(x, t) =m J f(x, v, t)dv, (8) 

pu(x,t) = m j f(x,v,t)vdv, 

f v 2 
pe(x,t)=m / f(x,v,t)—dv 

where m is the atomic/molecular mass. Supplementing 
these formal integrations with additional assumptions 
(e.g., small deviations from local thermodynamic equi- 
librium), one finally arrives at the desired equations for 
the macroscopic observables, typically the Navier Stokes 
equations of fluid dynamics (for the time being, we re- 
strict ourselves to the case of isothermal fluids, for which 
the energy equation is not needed): 

dtp + div(pu) = 0, (9) 
dtpu + div(puu) = —VP + div(2^Vtt + Aldivw), 

where p is the fluid density, u the fluid speed, P the 
fluid pressure, the overline denotes the symmetric dyad 
Vtt+(V«) T , the superscript T indicates the transpose, p, 
and A are the shear and bulk dynamic viscosities respec- 
tively, and 1 denotes the tensor identity. The Navier- 
Stokes equations keep no track of the discrete nature of 
the underlying microscopic world, and are paramount for 
the quantitative description of macroscopic systems. 



F. The top-down approach 

Each jump down the BBGKY ladder removes irrele- 
vant information so that, in the end, the 10 23 atomistic 
trajectories are replaced by the evolution of a handful of 
continuum hydrodynamic fields. 



The BBGKY approach is formally elegant and very 
fruitful for further theoretical insight and analysis. In 
fact, it is perfectly positioned to borrow the powerful 
mathematical machinery of classical and quantum sta- 
tistical field theory, such as perturbative methods, dia- 
grammatic techniques and the like. Less noted, perhaps, 
is the fact that the resulting equations prove exceedingly 
difficult to solve in actual practice. This is true even at 
the coarsest level: the Navier- Stokes equations are no- 
torious for posing one of the hardest problems left in 
classical (i.e., non quantum) physics, namely fluid turbu- 
lence. It makes sense therefore to think of complementary 
routes to BBGKY, more aligned with the spirit of model- 
building and computational tractability rather than with 
amenability to analytical treatment. 

An emerging and still fast-developing strategy along 
these lines is provided by fictitious dynamics methods. 
The idea is to introduce effective molecules ( "computa- 
tional quasi-particles"), each representing a huge num- 
ber, say R, of real ones, so that the number of effective 
molecules we must deal with is no longer of the order of 
Avogadro's number, but of order Nr = Na/R « Na 
instead. In terms of these effective molecules, the New- 
tonian equations are as follows: 

M I ^-=F I (X)+D I , I=1,N R (10) 

where Xj represents a coarse-grained coordinate, and 
Mj = m i i s the total mass of the effective "macro- 

molecule". The term Di collects all the details of the 
underlying fine-scales and disappears only in the trivial 
case of a linear dependence of the force Fj on the molec- 
ular coordinates. Realistic forces are generally inverse 
powers of the intermolecular distance and, consequently, 
Dj 7^ 0. A proper renormalization/closure procedure 
would attempt to incorporate the effects of the fine scales 
into appropriate (and most likely very complicated) ex- 
pressions for a renormalized force Fj = Ft + CDi, C 
being some form of projection operator making the renor- 
malized force available in terms of the coarse-grained co- 
ordinates. Most simply, one sets Dj = and proceeds 
with the solution of the same Newtonian equations, only 
applied to a much smaller set of molecules: this is the 
very successful path taken by Molecular Dynamics (Alder 
and Wainwright, 1959). 

However, once it is agreed that the ultimate aim is 
macroscopic physics, say solving the Navier-Stokes equa- 
tions, even Molecular Dynamics is still redundant of 
needless micro- details. We may want more than just set- 
ting Dj = 0, and look for instance for expressions of 
the coarse-grained force Fj that are simpler than those 
associated with true intermolecular interactions. The 
principle of Occam's razor, i.e., maximum simplicity im- 
plies that one should choose the simplest coarse-grained 
dynamics compatible with the target macroscopic equa- 
tions. This is the hard-core idea of "minimal molecu- 
lar dynamics" : relinquish as many microscopic details as 
possible right at the outset (atomistic level), making sure, 
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however, that the basic symmetries, conservation laws 
and evolutionary constraints needed to ensure the cor- 
rect macroscopic behaviour are preserved in the process. 
As to the first principle of thermodynamics, this means 
conserving all- and only-the microscopic invariants (mass, 
momentum, energy), while compliance with the second 
principle of thermodynamics imposes the existence of a 
suitable monotonically increasing function of time, the 
Boltzmann ii-function and the related entropy. The ac- 
tual realization of the top-down approach is by no means 
unique (for a recent form of dissipative particle dynam- 
ics, see Hoogerbrugge and Koclman, 1992, Espanol and 
Warren, 1995, Flekkoy and Coveney, 1999). In the fol- 
lowing, however, we shall refer to Lattice Gas Cellular 
Automata (LGCA), a particularly appealing instance of 
minimal particle dynamics developed in the mid 1980's, 
which provided the roots of the Lattice Boltzmann (LB) 
method. 



III. LATTICE GAS CELLULAR AUTOMATA 

The theory of lattice Gas Cellular Automata is a rich 
subject and has been recently made accessible in full de- 
tail by a scries of beautiful monographs (Rothman and 
Zaleski, 1997, Chopard and Droz, 1998, Rivet and Boon, 
2001). Here we shall take a substantial shortcut and pro- 
ceed by example. 

A. The Frisch-Hasslacher-Pomeau automaton 

Let us begin by considering a regular lattice with 
hexagonal symmetry such that each lattice site is sur- 
rounded by six neighbors identified by six connecting 
speeds Vi = v iai i — 1,6, the index a = 1,2 running 
over the spatial dimensions x, y (see Figure 2, which also 
includes a rest-particle with zero speed) . Each lattice site 
hosts up to six particles with the following prescriptions 

• All particles have the same mass m = 1 

• Particles can move only along one of the six direc- 
tions defined by the discrete displacements v ia . 

• In a time-cycle (set at unity for convenience) the 
particles hop synchronously to the nearest neigh- 
bor in the direction of the corresponding discrete 
vector Vi a . Longer or shorter jumps are both for- 
bidden, which means all lattice particles have the 
same energy. 

• No two particles sitting on the same site can move 
along the same direction v ia (Exclusion Principle). 

These prescriptions identify a very stylized gas analogue, 
whose dynamics is made purposely unaware of micro- 
scopic details of real-molecule Newtonian dynamics. In 
a real gas molecules move along any direction (isotropy) 
whereas here they are confined to a hexagonal cage. Also, 



real molecules can move virtually at any (subluminal) 
speed, whereas here only six mono-energetic beams are 
allowed. Amazingly, this apparently poor representation 
of true molecular dynamics has all it takes to simulate 
realistic hydrodynamics! With the prescription kit given 
above, the state of the system at each lattice site is un- 
ambiguously specified in terms of a plain yes/no option 
saying whether or not a particle sits on the given site. 
This dichotomic situation is readily coded with a single 
binary-digit (bit) per site and direction so that the entire 
state of the lattice gas is specified by 6N bits, N being 
the number of lattice sites. Borrowing the parlance of 
second quantization, we introduce an occupation number 
rii, such that 

«<(*,*) = {1,0} (11) 

depending on whether or not the lattice site x hosts a 
particle with speed Vi at time t. The collection of oc- 
cupation numbers rii(x,t) over the entire lattice defines 
a 6N dimensional time-dependent boolean field whose 
evolution takes place in a boolean phase-space consist- 
ing of 2 6N discrete states. This boolean field takes the 
intriguing name of Cellular Automaton to emphasize the 
idea that not only space and time, but also the dependent 
variables (matter) take on discrete (boolean) values. The 
fine-grain microdynamics of this boolean field cannot be 
expected to reproduce the true molecular dynamics to 
any reasonable degree of microscopic accuracy. However, 
as is known since Gibbs, many different microscopic sys- 
tems can give rise to the same macroscopic dynamics, and 
it can therefore be hoped that the macroscopic dynamics 
of the lattice boolean field would replicate real-life hy- 
drodynamic motion even if its microdynamics does not. 

B. Boolean microdynamics 

Let us now prescribe the evolution rules of our cellular 
automaton. Since we aim at hydrodynamics, we should 
address two basic mechanisms: 

- Free- streaming 

- Collisions 

Free streaming consists of simple particle transfers 
from site to site with discrete speeds V{. Thus, a par- 
ticle sitting at site x at time t with speed Vi a will move 
to site x + Vi at time t + 1. 

In equations: 

rii(x + v it t + 1) = rii(x,t) (12) 

This defines the discrete free-streaming operator Si as 

Sini = ni(x + Vi,t+l)-ni(x,t) (13) 

This operator is a direct transcription of the Boltzmann 
free-streaming operator, D t = d t + v a d a , to a discrete 
lattice in which space-time variables are discretized ac- 
cording to the synchronous "light-cone" rule: 
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Ax h 



Via At. 



(14) 



Once on the same site, particles interact and reshuffle 
their momenta so as to exchange mass and momentum 
among the different directions allowed by the lattice (see 
Figure 3). 

This mimics the real-life collisions taking place in a real 
gas, with the crude restriction that all pre- and post- 
collisional momenta are forced to "live" on the lattice. 
As compared with continuum kinetic theory, Lattice Gas 
Cellular Automata introduces a very radical cut of de- 
grees of freedom in momentum space: just one speed 
(all discrete speeds share the same magnitude c = 1, 
hence the same energy) and only six different propaga- 
tion angles. Not bad for an original set of six-fold infi- 
nite degrees of freedom! Space-time is also discretized 
(see (p.4[)), but this is common to all computer simu- 
lations of dynamical systems. At this stage, it is still 
hard to believe that such a stylized system can display 
all of the complexities of fluid phenomena. And yet it 
does! The reader acquainted with modern statistical 
mechanics smells the sweet scent of universality: for all 
its simplicity, the Frisch-Hasslacher-Pomeau automaton 
may display the same large-scale properties of a real fluid 
(Kadanoff, 1986), such as propagation of sound waves, 
vortex interactions and energy dissipation. The name of 
this magic is symmetry and conservation. 

Let us dig a bit deeper into this matter. 

Let us consider the Frisch-Hasslacher-Pomeau collision 
depicted in Figure 3: Albeit stylized, this collision shares 
two crucial features with a real molecular collision: 

• It conserves particle number (2 before, 2 after) 

• It conserves total momentum (0 before, after) 

Symbolically, its effect on the occupation numbers is a 
change from n.$ to n[ on the same site 



Ci{n) 



(15) 



where n = [m, Tii . . . Tib] denotes the set of occupation 
numbers at a given lattice site. 

To sum up, the final Lattice Gas Cellular Automata 
update rule reads as follows: 



SiTli — Ci 



(16) 



or, which is the same: 



m(x + Vi,t + 1) = n'iixjt) 



(17) 

where all quantities have been defined previously. The 
equations (16[ |l7|) represent the microdynamic equation 
for the boolean lattice gas, the analogue of Newton equa- 
tions for real molecules. This equation constitutes the 
starting point of a Lattice BBGKY hierarchy, ending up 
with the Navier-Stokes equations. At each level, one for- 
mulates a lattice counterpart of the various approxima- 
tions pertaining to the four levels of the hierarchy (see 
Figure 1). 



C. Merits and pitfalls of Lattice Gas Automata 

The major appeals of Lattice Gas Cellular Automata 
are: 

• Round- off free computing (boolean algebra is exact) 

• Memory savings ( only one bit per degree of free- 
dom) 

• Virtually unlimited potential for parallel computing 

These points fueled the excitment around Lattice Gas 
Cellular Automata as a potentially revolutionary tool for 
computational fluid dynamics (and more). In particu- 
lar, the hope was to attack the infamous problem of fluid 
turbulence. Let us just mention that this relates to the 
dynamics of flows where dissipative effects are very small 
as compared with advection. The relative strength of ad- 
vection versus dissipation is measured by a dimensionless 
number known as the Reynolds number: 



Re = 



UL 



(18) 



where U is a typical flow speed on a macroscopic scale L 
(the size of the device) and v is the kinematic viscosity 
of the fluid. Based on dimensional scaling theories (A.N. 
Kolmogorov, 1941), it can be shown that the number of 
degrees of freedom associated with a turbulent flow at a 
given Reynolds number Re is approximately i?e 9 / 4 . Since 
Re ~ 10 6 is a commonplace in daily life (that's more or 
less what we experience by driving our car at cruising 
speed of about 100 Km/h in full compliance with traffic 
regulations...), the simulation of such flows implies the 
solution of about 10 14 degrees of freedom: less than Avo- 
gadro's number, to be sure, but still too much for any 
foreseeable computer. These figures say it all as to the 
hunger for innovative mathematical and numerical meth- 
ods for fluid turbulence! Unfortunately, closer inspection 
into the details of the Lattice Gas Cellular Automata 
method reveal a number of difficult problems: 

• Statistical noise 

• Complexity of the collision operator 

• Small number of collisions 

Statistical noise relates to the fact that in order to 
extract a smooth hydrodynamic signal, averages over 
many boolean variables are required, thereby eclipsing 
the memory savings provided by the boolean micrody- 
namic representation. 

Exponential complexity relates to the exponential es- 
calation of the collision operator as more physics is added 
to the model, say, more than one fluid species, or simply 
by moving to higher dimensions. The problem is that the 
complexity of the corresponding collision operator grows 
exponentially, roughly as 2 b , where b is the number of bits 
per site, so that the original simplicity is rapidly lost. 
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Low collisionality is also related to the paucity of dis- 
crete speeds. Only a relatively small fraction of phase 
space is collisionally active since many collisions are sim- 
ply not compatible with conservation principles. Few 
collisions mean long mean-free-path, hence high momen- 
tum diffusivity v 1 hence low-Reynolds numbers, and the 
dream of simulating highly turbulent flows fades away. 

In spite of the remarkable progresses achieved in the 
late 1980's, the Lattice Gas Cellular Automata algo- 
rithms remained rather heavy and stiff with respect to 
the collision rules. On the other hand, these inconve- 
niences were not compensated by a dramatic advantage 
in terms of accessible Reynolds numbers. As a result, the 
interest in Lattice Gas Cellular Automata as a tool for 
high Reynolds flow simulations leveled off in the early 
1990's. This was precisely the time the Lattice Boltz- 
mann method took off. 



IV. LATTICE BOLTZMANN EQUATIONS 

The theory of the lattice Boltzmann equation (LBE) 
(for a review see, Benzi, Succi, Vegassola, 1992, Qian, 
Succi, Orszag, 1995 and Chen, Doolen, 1998) begins with 
a straightforward floating point recast of the boolean evo- 
lution equation of Lattice Gas Cellular Automata dynam- 
ics (Frisch, Hasslacher and Pomeau, 1986; Frisch et ai, 
1987). With reference to a set of b populations, fi(x,t), 
with {?' = 1, . . . , b} representing the probability for a par- 
ticle to reside on a lattice site x at time t with discrete 
velocity Vi, the LBE reads as follows: 

fi(x + v t ,t + 1) = fi(x, t) + Ci[h, ...,/(,], (19) 

where Ci[f] is the collision operator, a polynomial of de- 
gree b constructed explicitly out of all allowable n-body 
collisions, 2 < n < b, among populations sitting on the 
same site x at time t. The above LBE has a simple 
interpretation: Particle population fi(x + m, t + 1) is 
equal to the post-collision f[(x,t) value advected from 
the "upper wind" position x at the previous time t. Here 
fl = fi + Ci. Thus Ci represents the change of the parti- 
cle population by collisions. Compliance with mass and 
momentum conservation laws imposes the following con- 
straints on the collision operator: 

i i 

If there is an energy degree of freedom, then we also have: 
5>C<=0, (21) 

i 

where Ci = vf /2 for an ideal gas. As in continuum kinetic 
theory, this collision operator admits a detailed balance 
condition when all populations of different particle veloc- 
ities are in equilibrium, {/,; = /? q , Vi}: 

a[/ eq ]=0 (22) 



Regardless of the detailed collision processes, the solution 
of Eq. (|22| ) takes on a generic form dictated by the basic 
conservation laws, 

fP = ^ (23) 

where li = A + B • Vi + Cti is a linear combination of col- 
lisional invariants, while A and B, in turn, arc functions 
of local hydrodynamic quantities, 

p(x,t) = ^2f i (x,t) (24) 



pu(x,t) =^2vifi(x,t) (25) 

i 

Once again, for systems with an energy degree of free- 
dom, the total energy can be defined as, 

p(DT + u 2 )/2 = J2^ifi (26) 

i 

where D is the number of spatial dimensions and T can 
be interpreted as the temperature of the fluid. 

The LB equation ( |l9| ) obeys the so-called semi-detailed 
balance property for the collision operator as its Lat- 
tice Gas Cellular Automata predecessor (Frisch, Hass- 
lacher and Pomeau, 1986; Frisch et ai, 1987). It can be 
shown that such a kind of collision operator admits a lo- 
cal H theorem with the following discrete Boltzmann en- 
tropy function, h-Q(x,t) = '^2 i fi(x,i)ln fi(x,t). In other 
words, an H theorem can be defined on each local lattice 
site, so that, without external disturbances or boundary 
influences, /ib is a non-increasing function of t satisfying 
the dynamics of Ci. Furthermore, the minimum value of 
liB(x,t) is attained when the particle populations {/$} 
assume the equilibrium form given by (^3|), where A, B 
and C are constants determined by the values of mass, 
momentum and energy. Borrowing from the terminology 
of nonlinear dynamic systems, the existence of an H the- 
orem ensures that the equilibrium distribution is not only 
a fixed point solution of ( j2^ ) but also an attractor of the 
collisional dynamics. On the other hand, this is a weaker 
H theorem compared to its continuum counterpart, for 
the above does not include advection of particles among 
lattice sites according to (|l^) . We will realize in a subse- 
quent section that, unlike in continuum Boltzmann dy- 
namics, the local H theorem in Lattice Gas Cellular Au- 
tomata given above does not (except for the isothermal 
case, in which energy conservation is replaced by a con- 
stant temperature-like parameter) automatically lead to 
a global one, referred to as the global i7-theorem (Frisch, 
Hasslacher and Pomeau, 1986; Frisch et ai, 1987, Chen, 
1995; Chen, 1997). Generally the global H function is 
given by, 

H(t) =J2h B {x,t). 
x 
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Existence of a global H theorem in LB guarantees an 
asymptotically homogeneous spatial distribution of par- 
ticles as time t — ► oo, and hence it provides a well-defined 
global hydrodynamic stability. As in standard kinetic 
theory, hydrodynamics concerns properties around the 
local equilibrium, 

fi = /r + /r (27) 

where the non-equilibrium component, /" G , is supposed 
to scale like kf* q , where the smallness parameter k, 
known as the Knudsen number, is the ratio of the parti- 
cle mean free path to a typical hydrodynamic scale. Or- 
dinary fluids feature Knudsen numbers below 0.01, see 
Chapman and Cowling, 1970.) The existence of a lo- 
cal H theorem makes a perturbation analysis purposeful 
since the dynamic system always approaches a defined 
equilibrium / 8 eq and evolves in the neighborhood of such 
a defined point. Perturbative treatment of the LBE up to 
second order in k (Chapman-Enskog procedure) is then 
expected to yield the hydrodynamic equations 

dtp + d a pu a = 
dtpu a + d b P ab = (28) 

where 

p ab = ^> ia <wr +5> fa <wr (29) 

i i 

= gpu a u b + PS ab + pv\d a u b + d b u a - -^{d c u c )5 ab \ 

is the momentum-flux tensor comprising both non- 
viscous and dissipative components, the latter being pro- 
portional to the fluid viscosity v. In the above, P is the 
resulting pressure, and D is the dimensionality of a lat- 
tice. Note that the prefactor g in the advection term 
which signals a potential breaking of Galilean invariance 
whenever j ^ 1. On the other hand, since g depends 
only on the fluid density p, for incompressible flows in 
which the density is constant in space and time, Galilean 
invariance can be recovered by a simple rescaling of time, 
pressure and viscosity: 

t' = gt, i/ = u/g, P' = P/g 

Clearly, even within the framework of incompressible 
flows, this rescaling does not extend to more general sit- 
uations, such as multi-component or multi-phase flows. 

A. The fully non-linear LBE 

The fully non- linear LBE, cq. (JTJ), was first pro- 
posed by McNamara and Zanetti (McNamara, Zanetti, 
1988), who realized its potential for doing away with 
the statistical noise problem affecting Lattice Gas Cellu- 
lar Automata simulations (Orszag, Yakhot, 1986, Succi, 
Santangelo, Benzi, 1989). All other weaknesses still re- 
mained. In general, the factor g is not equal to unity, 



indicating a violation of Galilean invariance. Though a 
choice of sufficient lattice symmetry ensures a rotation- 
ally invariant form of P ab , the resulting local equilibrium 
in typical Lattice Gas Cellular Automata collisions does 
not have a suitable functional form for achieving a cor- 
rect Navier-Stokes hydrodynamics, which requires g = 1. 
Indeed, in order to obtain the correct form of the Navier- 
Stokes equation, the local equilibria must comply with 
a specific form as a function of the local hydrodynamic 
variables, p and u, (S. Chen et ai, 1991; H. Chen, S. 
Chen and Matthaeus, 1992; Qian, d'Humieres and Lalle- 
mand, 1992). 

In addition to incorrect hydrodynamics, Lattice Gas 
collisions have problems in attaining low viscosity. The 
fluid viscosity is approximately given by v = 3DIvt, 
where I is the particle mean-free path and vt the typical 
particle thermal speed. In many problems, for instance 
turbulence, we are interested in very low-viscous flows, 
which means very short mean-free paths. On the other 
hand, short mean-free paths imply many collisions. 

We are thus forced to consider as far as possible all 
allowable types of collisions among the set of discrete ve- 
locities. Unfortunately, it turns out that the full collision 
operator Ci involves an exponential barrier of the order 
of 2 b operations, in direct contradiction with the basic 
commitment to simplicity and computational efficiency 
of the method. 



B. The LBE in scattering form 

Building upon the idea that many-body collisions are 
not essential to achieve the correct hydrodynamic limit, 
Higuera and Jimenez (Higuera and Jimenez, 1989) re- 
alized that the collision operator can be reduced to a 
dramatically simpler two-body scattering expression: 

c -X><'/, ( 3 °) 

3 

where the scattering matrix Ay is basically the Ja- 
cobian of the fully non-linear operator Cj evaluated 
at the uniform equilibrium values /, = p/b. The 
above expression turns a daunting 2 b complexity into a 
much more manageable b 2 one, thus opening the way 
to three-dimensional Lattice Boltzmann hydrodynamics. 
Of course, compliance with an _£/-theorem is no longer 
guaranteed because the local equilibrium is no longer the 
direct result of collisional dynamics. 

In addition, since the scattering matrix Ay is still re- 
lated one-to-one to the underlying Lattice Gas Cellu- 
lar Automata microdynamics, the corresponding Lattice 
Boltzmann equation shares the same limitations in terms 
of high viscosity, i.e., low Reynolds numbers. 
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C. The self-consistent LBE 

This last limitation can be wiped out by a mere change 
in perspective. Instead of deriving the Navier-Stokcs 
equation bottom-up (here bottom means the atomic 
level) from a truly TV-body discrete dynamical system, 
we can construct it top-down from the sole requirement 
of compliance with the Navier-Stokes equations (Higuera, 
Succi and Benzi, 1989, Succi, Bcnzi, Higuera, 1991). The 
idea is to recognize that, as far as hydrodynamics is 
concerned, the key notions of scattering matrix and lo- 
cal equilibrium can be prescribed at the outset instead 
of being derived from an underlying (discrete) micro- 
dynamics. Mathematically, this amounts to prescribing 
the scattering matrix in spectral-form: 

^• = 5> fc i#\ (31) 
k=i 

where projects along the fc-th eigenvector in kinetic 
space and Aj. is the corresponding eigenvalue (A& = 
for conserved quantities). Of particular interest is the 
leading non-zero eigenvalue, which is in direct control of 
the fluid viscosity v <~ A -1 . 

Once this point of view is endorsed, namely that lo- 
cal equilibria can be "freely" chosen within the conserva- 
tion constraints, and that the scattering matrix can be 
selected a priori on the sole basis of conservation con- 
straints, nothing prevents this Lattice Boltzmann equa- 
tion from attaining as low viscosity as the lattice dis- 
creteness permits. Full contact with classical computa- 
tional fluid dynamics is made (Chen H., S. Chen, W. 
Matthaeus, 1992). More importantly, the top-down ap- 
proach, which proved so fruitful and influential for sub- 
sequent developments of LB theory, is established. 

D. The lattice Bhatnagar Gross-Krook equation 

The LBE story has still one more nugget in store. 

The scattering matrix can also be viewed as a mul- 
tiple scale relaxation operator, one scale for each non- 
zero eigenvalue. Since we are basically interested in a 
single transport parameter, the fluid viscosity, a single 
eigenvalue should do. Indeed we can replace the full ma- 
trix Aij with a single-parameter diagonal form —uiSij, 
describing a single-time relaxation around a prescribed 
local equilibrium /° q . In its simplest and by now most 
popular form, the relaxation-approximation corresponds 
to the following Lattice BGK (Bhatnagar-Gross-Krook) 
equation (Bhatnagar, Gross and Krook, 1954): 

fi(x + Vi,t+l)-fi(x,t) = -uj[fi(x,t)-fP(x,t)], (32) 

where exact Navier-Stokes hydrodynamics is obtained 
with the following choice of local-equilibrium form, 

/r=^p|i+^+ - 2T2 — y (33) 



where the symbol : stands for tensor scalar product. In 
the above, is a suitable lattice-dependent weighting 
factor, and the temperature T = 1/3 in typical isother- 
mal lattice BGK models (S. Chen et at, 1991; H. Chen, S. 
Chen and Matthaeus, 1992; Qian, d'Humicrcs and Lalle- 
mand, 1992; Chen, Teixeira and Molvig, 1997). Since 
temperature is frozen to a constant value, these lattices 
models should best be denoted as "a-thermal" rather 
than "iso-thermal" . More comments on this delicate is- 
sue shall be presented later in this paper. For a thor- 
ough and beautiful discussion of the subtleties of lattice 
thermo-hydrodynamics, the reader is recommended to 
consult Rivet & Boon, Chapter IV. 

It is readily recognized that this formulation leads to a 
Galilean invariant Navier-Stokes equation (up to terms of 
order Mac/i 4 , where Mach is the Mach number, namely 
the ratio of fluid to sound speed), for a fluid of viscosity 
v <~ T/uj. The LBE scheme in relaxation form has met 
with significant success in the last decade in simulating 
a variety of fluid flows. Indeed, the lattice BGK equa- 
tion and subsequent straightforward extensions to allow 
a variable Prandtl number (ratio of momentum to heat 
diffusivity) (Keizer, 1987, Chen et al 1997), represents 
the method of choice in the field. Very interesting vari- 
ants, combining the best features of the Lattice Boltz- 
mann equation in scattering form and Lattice BGK have 
also been developed (D'Humicrcs, 1992). 

A few years later it was shown that the Lattice BGK 
can be derived from the continuum Boltzmann BGK 
equation by a Grad-like moment expansion supplemented 
with numerical quadrature for the actual evaluation of 
the kinetic moments (He and Luo, 1997; Shan and He, 
1998). Realizing such a connection might prove very 
useful for establishing new lattice BGK's starting from 
model continuum Boltzmann equations, hopefully includ- 
ing physics beyond the hydrodynamic level. 

While focus on the conservation laws (hydrodynamic 
constraints) was the leit motif of all the aforementioned 
developments, compliance with the second thermody- 
namic principle, namely the existence of the H theorem, 
was somehow left in a subdued light by the top-down 
approach. That neglect had no serious consequences, be- 
cause, thanks to a gracious smile from Lady Luck, being 
cavalier about an underlying H theorem proved very for- 
giving for isothermal LBE flows (the bulk of early LBE 
research) even at very low viscosities, where numerical 
stability is severely probed. By now, we have learned 
that this favorable behavior is due to the existence of 
an underlying H theorem, as we shall detail in the next 
section. 



V. H THEOREM IN DISCRETE PHASE-SPACE 

The H theorem is a milestone of non-equilibrium sta- 
tistical mechanics, since it provides a conceptual link be- 
tween the reversible laws of the microworld and the one- 
sided nature of macroscopic phenomena (Lebowitz, 1993; 
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Lieb, 2000). It is also a fundamental concept in com- 
putational physics, where compliance with an H theo- 
rem is often perceived as a byword for numerical stabil- 
ity (Perthame and Tadmor, 1991; Natalini, 2000; Junk 
and Klar, 2000). ff-compliant Discrete- Velocity Boltz- 
mann (DVB) equations have been known for a long time 
(Broadwell, 1964, Gatignol, 1965), but, as discussed pre- 
viously, they are not meant to compete with Navier- 
Stokes solvers because their main focus is maximum like- 
lyhood to the Boltzmann equation in the high-Knudscn 
regime. This leads to discrete collision operators which 
are just too complicated to serve as a practical tool for 
purely hydrodynamic purposes. 

Computational complexity dissolves by using relax- 
ation approximations, notably lattice BGK, and it is 
therefore imperative to explore the possibility of estab- 
lishing H theorems for lattice BGK. Before plunging into 
this theme, a few reflections on the classical case are in 
order. 



A. Reflections on the continuum case 

The classical Boltzmann H theorem is so familiar that 
it is worthwhile to look at it from the perspective of mod- 
eling, in order to better understand what is actually lost 
and needs to be rearranged when going from the contin- 
uum to the Lattice Boltzmann case. 

First, we find that the familiar local Maxwell distri- 
bution function (the "Maxwellian" ) minimizes the H 
function, H = J f In fdv, as soon as the five hydro- 
dynamic invariants (mass, momentum and energy), are 
fixed. Furthermore, we observe that all higher-order mo- 
ments of the local Maxwellian distribution are just in 
the right form for the Chapman-Enskog expansion to de- 
liver the Navier-Stokes equations. Backed with this in- 
formation, we start looking for a kinetic model acting 
as a continuous-time constrained minimization process 
of this H, in such a way as to conserve the hydrody- 
namic moments. We readily find several such models: 
the Boltzmann equation itself, the BGK model, the dif- 
fusion equation, and so on. All of these models differ 
only in the way relaxation to the local equilibrium takes 
place; some of them do require explicit knowledge of the 
local Maxwellian (the BGK), some others don't. They all 
deliver the Navier-Stokes equation in the end, with the 
only differences among them occurring in the transport 
coefficients, whose evaluation is simple for the BGK and 
requires some (non-negligible) work for the Boltzmann 
equation. In all models the collision integral has the lo- 
cal Maxwellian as its zero point. The rate at which the 
H function decreases in time is just equal to the entropy 
production (we recall that entropy is the space integral 
of the —H function) and entropy production becomes 
zero in the local Maxwell state. The local Maxwellian 
distribution is then characterized in three different but 
equivalent ways: It is the minimum of H , it is the zero 
point of the collision integral, and it is the zero point of 



the entropy production. 

When moving to the lattice world the basic question 
to be addressed is: how does the H theorem transform in 
the discrete-time easel 

The first good news is that all velocities are in some 
finite range, so that we can think of local equilibria for 
the given density and velocity alone, without necessar- 
ily including energy. In the (non-relativistic) continuum, 
energy must be included, for otherwise nothing would 
prevent integrals over velocities from diverging. The flip 
side is that since discrete speeds come in a very(!) fi- 
nite number, the Boltzmann H function does not work: 
For any known lattice, a straightforward computation 
demonstrates that the local Maxwellian equilibrium does 
not imply correct expressions for the higher-order mo- 
ments. On the other hand, a priori, there is a huge class 
of convex functions at our disposal, and in order to illus- 
trate the idea of what local equilibria may look like in the 
lattice context, we shall present an example of a solvable 
lattice entropy, i?i/a = 2i=i fiVTi 1 - F° r this entropy, 
the local equilibrium can be found explicitly (Karlin et 
al, 1998) 



2cf (l + Vl - M 2 ) 
(34) 

Here c 2 is the speed of sound squared, and M 2 = u 2 /c 2 is 
the local Mach number squared. The above equilibrium 
is Galilean-invariant only in the limit of vanishing flow, 
M — > 0. At any finite flow speed, a quadratic (in the 
Mach number) anomaly is apparent. Because advection 
terms in the Navier-Stokes equations are also quadratic 
in velocity, such a rapid growth of anomalies rules out 
the use of the function H\ji for constructing the Lat- 
tice Boltzmann method. This is the typical situation of 
first-generation Lattice Gas and Lattice Boltzman mod- 
els. Nevertheless, it is instructive to look at Eq. (|34|) and 
see how lattice equilibria differ from the local Maxwellian. 
Recall that the local Maxwellian is a well-defined func- 
tion for all values of the average velocity. In contrast, the 
local equilibrium ([54j) is positive only for M < 1, and it 
does not exist as a real- valued function, for M > 1 . This 
means that no collision mechanism is able to equilibrate 
the non-equilibrium deviations produced by the super- 
sonic motion, and therefore no macroscopic dynamics for 
this regime can exist. This conveys an intriguing flavour 
of relativistic mechanics, which is after all not surprising 
since lattice molecules move at the lattice speed of light 
c equal to the maximal length of the lattice link. 



1 This entropy function belongs to the class of Renyi (or Tsallis) 
entropies, H p = (fP +1 - f)/p, with p = 1/2. This is a glorious 
class of functions of modern statistical physics, for it provides 
the stepping stone for the celebrated replica method in spin-glass 
theory, multifractal ideas and also non-extensive thermodynam- 
ics (Tsallis, 1988). 
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In order to proceed sensibly, three basic requirements 
must be faced: 

• Galilean invariance 

• Readability 

• Solvability 

each of which we now discuss in some detail. 



B. Galilean invariance 

Galilean invariance requires kinetic equilibria to de- 
pend on the relative speed v— u ("peculiar" speed) rather 
than on the absolute molecular speed v itself. In the con- 
tinuum, this is ensured by the Max wcllian d ependence, 
~ exp{ — (v— u) 2 /vt}, where vt = y^fceT/m is the ther- 
mal speed, which sets the natural scale for molecular fluc- 
tuations around the fluid speed u. In the Lattice Boltz- 
mann setting, we have already committed ourselves to 
low Mach numbers and no compressibility effects. There- 
fore, it seems reasonable to proceed with expansion of the 
local Maxwellian around the global equilibrium (u = 0). 
Since the Maxwellian is a transcendental function, large 
departures from global equilibrium require virtually an 
infinite number of terms of this expansion. 

Terms that correspond to kinetic excitations on top of 
the uniform "ground state" are described by higher or- 
der polynomials in the velocity variable. Since a finite 
set of discrete speeds can only support a finite number 
of these excitations, breaking of Galilean invariance can- 
not be avoided. The above considerations suggest that 
Galilean invariance can be recouped to some order if lo- 
cal equilibria can be expressed in the form of finite-order 
polynomials of the flow speed. 

The question is then whether suitable polynomials can 
be found that still comply with a discrete H theorem. 
The answer to this question appears to be negative (Wag- 
ner, 1998): There are no convex entropy functions whose 
local equilibria are the polynomials of the form used in 
the lattice BGK. 

In analogy with lattice field theory, we shall call 
Galilean-invariant discrete entropies as perfect, in that 
they "hide" the underlying lattice discreteness. Are 
there perfect entropy functions for the Lattice Boltzmann 
method? At the time of this writing, no perfect lattice 
entropy is known, and most probably there is none, so 
that we are led to conclude that the Lattice Boltzmann 
method is not capable of reproducing the full properties 
of the continuum Boltzmann equation. 

However, quasi-perfect entropies, namely entropies 
which are not affected by lattice discreteness (up to 
fourth order terms in the Mach number) have indeed 
been found by a customized, lattice-dependent procedure 
(Karlin, Ferrante and Ottinger, 1999). As an illustration, 
for a three-state one-dimensional lattice with Vq = and 



v T i = TL the following quasi-perfect Boltzmann-like H- 
function is identified: 

H = f ln(/ /4) + /+ In /+ + /_ In /_ . (35) 

Given the fact that the Lattice Boltzmann equation is 
itself a second order approximation in the Mach num- 
ber to the Navier-Stokes equations, these quasi-perfect 
entropies must be regarded as definitely adequate for hy- 
drodynamic purposes. 

The explicit form of the local equilibria corresponding 
to the quasi-perfect Boltzmann-like entropy functions is 
not known, in general. However, polynomial approxima- 
tions can be found up to the relevant order in Mach num- 
ber. These approximations coincide with those estab- 
lished earlier for the lattice BGK, and this now explains 
why the Lattice Boltzmann works in this case: These spe- 
cific polynomial equilibria survived among other possible 
Galilean-compliant polynomials because they are com- 
putationally more stable than others, and it is precisely 
these polynomials which are supported by the quasi- 
perfect entropy functions. The existence of quasi-perfect 
lattice entropies achieves a formal compliance with the 
second law of thermodynamics and provides a strong 
incentive to pursue further development of the Lattice 
Boltzmann method based on the entropy maximization 
principle. 

C. Readability 

Realizability is simply the condition that local equilib- 
ria resulting from an entropy maximization procedure be 
real-valued and between zero and one: 

< /r < 1 (36) 

By itself, realizability is neither a necessary nor a suf- 
ficient condition for stability, although it generally helps 
stability (Rcnda et al, 1998). At any rate, it is a good 
pre-requisite, at least until we learn to deal sensibly with 
negative distribution functions. 

It should be appreciated that gradient-driven depar- 
tures from local equilibrium can violate the realizability 
constraints even if the local equilibria don't. Therefore 
the stability domain in kinetic space is by necessity a 
subset of the realizability domain. Manifestly, the poly- 
nomial nature of the discrete local equilibria is a poten- 
tial danger to the realizability constraint, a danger which 
simply does not exist in the continuum. 

Description of realizability domains on the lattice is 
a nontrivial problem with an interest of its own. In a 
recent work, Boghosian's et. al (2001) applied the pow- 
erful Fourier- Motzkin technique of linear programming to 
classify automatically the realizability domain associated 
with a given set of constraints. 

This may help the systematic search for optimal lattice 
entropies. 
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D. Solvability 

Solvability refers to the possibility of expressing local 
equilibria as explicit functions of the conserved hydrody- 
namic variables, i.e., density p and flow speed u. This is 
very important since it permits one to explicitly recast 
kinetic theory in terms of an equivalent set of partial dif- 
ferential equations for space-time dependent continuum 
fields. Solvability has also a considerable practical impact 
on the efficiency of lattice BGK schemes, for it implies 
that local equilibria can be encoded once and for all as 
analytic functions of the hydrodynamic variables. Local 
equilibria resulting from non-solvable entropies must be 
recomputed numerically at each time step by iterative 
procedures and deprive lattice BGK of (part of) its sim- 
plicity and efficiency. A way out of this problem is to 
restore the Boltzmann-like collision operators which can 
be constructed from just the knowledge of the entropy. 
(We recall that even the single relaxation time property 
is not exclusively the property of the Bhatnagar-Gross- 
Krook equation) . This has been done recently ( Ansumali 
and Karlin, 2000). 



E. Discrete-time effects and the mirage of zero-viscosity 

Perhaps the most interesting feature of the discrete- 
time kinetics is that the if-theorem is no longer the 
same (qualitatively) as in the continuous-time case. The 
way the iJ-theorem features itself in the Discrete Veloc- 
ity Models, where time is still continous, is basically the 
same as in the classical theory. When time is discrete, 
things change considerably. 

The major distinction is that, in order to achieve in- 
terestingly low values of the transport coefficients, the 
lattice relaxation dynamics must proceed in artificially 
long (hence more effective) jumps going across the max- 
imum entropy point (equilibrium) , in a sort of two-sided 
over-relaxation process, which bears little resemblance to 
the smooth relaxational trajectory of the continuum case. 
The problem is best illustrated in geometrical terms (see 
Figure 4): If / is the set of populations at time t, then 
the collision operator gives the direction, C(f), in the 
kinetic space in which the populations must be moved 
(these directions change at different lattice nodes). Let 
us consider populations / + f3C[f], where (3 > 0. When 
tracing the ray / + j3C[f], starting with (3 = 0, we see 
that the function H is decreasing (because the classical 
entropy production, — ^\ Ci[f]dH/dfi, is positive), then 
comes to the minimum at some /?', then starts increasing 
which would not be allowed if time were continous. In 
discrete time we can still safely place the populations in 
the increasing branch until H (/ + 0C[f]) < H(f)), and 
finally come to a value (3* where the entropy just equals 
the initial value H(f). Formally, the condition: 

H(f)=H(f + f3*C[f}), (37) 



sets the limit for the over-relaxation. It has been shown 
(Karlin et al, 1998; Karlin, Ferrante and Ottinger, 1999; 
Boghosian et al, 2001) that this estimate reduces to the 
so-called linear stability interval, < uj < 2 for the lattice 
BGK, required by the positivity of transport coefficients, 
as soon as the state is close to the local equilibrium. Re- 
call that in the classical continuous-time kinetic theory, 
positivity of transport coefficients is ultimatively related 
to positivity of the local entropy production. It is re- 
markable that the lattice local H -theorem establishes the 
same result for the discrete-time case since positivity of 
transport coefficients follows now from the estimate 
In other words, it is clear that in the over-relaxation 
scenario there is no guarantee that the post-collisional 
state generally attains a smaller H-vahie than the pre- 
collisional state. This is controlled by the global shape 
of the entropy function far from equilibrium. In fact, if 
the value of u> is too large, the end-point may run into 
a lower entropy state, or some populations may even be- 
come negative. Both such outcomes set the stage for the 
type of kinetic instabilities hindering the application of 
the Lattice Boltzmann method to very low viscous, high 
Reynolds flows. 

This discussion shows that the "mirage" of zero- 
viscosity, so crucial for turbulence studies, falls within 
the general problem of finding a systematic formulation 
of the dynamics of dissipative systems far from equilib- 
rium (Prigogine, 1962; Ruelle, 1999). 

With the quasi-perfect entropy at hand, we make full 
contact with the second law of thermodynamics, and 
establish the "smallest" Boltzmann system: It has the 
properly defined local equilibrium, the H theorem, and 
the correct hydrodynamics to the minimal required order 
of approximation. 

The discrete-time f/-theorem suggests the possibility 
of increasing stability at low viscosity by implementing 
the estimate (|37|). In fact, Lattice Boltzmann schemes en- 
dowed with quasi-perfect entropies, and with the entropy 
estimate ( |37| ) implemented, are generally found to exhibit 
better numerical stability. Typical instabilities associ- 
ated with the lack of the _ff-theorem are shown in Figures 
5 and 6. In Figure 5 we show the density profile for a 
one-dimensional front in a shock tube at time t = 500 
(in lattice units) for the Lattice Boltzmann with entropy 
function (Ansumali and Karlin, 2000) (top picture, la- 
belled ELBM), the Lattice BGK model based on the 
polynomial ansatz of Qian, d'Humieres and Lallemand 
(1992), (middle picture labelled LBGK), and the Lattice 
Boltzmann model by Qian, D'Humieres and Lallemand, 
1991, (bottom figure, labelled LBE). The viscosity is set 
to v = 10 _1 /3. From Figure 5, we notice that, although 
LBE is producing the most accurate solution almost ev- 
erywhere, only ELBM is free of non-physical oscillations 
(known as "Gibbs phenomenon" in numerical analysis), 
typical of non-entropic numerical schemes and often a 
precursor of numerical instabilities. Indeed, a minor de- 
crease of the viscosity is found to disrupt the stability of 
the LBE simulation. By further decreasing the viscosity, 
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LBGK also becomes unstable, whereas ELBM does not. 
On the other hand, the better ELBM stability comes at 
the expense of some oversmoothing of the fronts, as it 
is visible by comparison with the thin line (top figure) 
giving the exact solution. The Figure 6, which shows 
the corresponding velocity profile, tells essentially the 
same story: the long-standing conflict between stability 
and numerical diffusion, a sort of "numerical uncertainity 
principle" (Boris, 1989). For very recent improvements 
of the ELBM approach in the direction of decoupling the 
conflicting issues of numerical stability and accuracy, see 
Ansumali and Karlin, (2002a, 2002b). 

F. Over-constrained equilibria 

We emphasize that local equilibria of the quasi-perfect 
entropy are "classical" : They minimize the entropy un- 
der constraints provided by the locally conserved fields, 
whereas the desired conditions for the non-conserved 
fields (the local equilibrium momentum flux tensor in the 
above consideration) follow from this solution. 

What happens if the set of constraints is enlarged in 
such a way as to include higher-order moments? 

It should be noted that this is a rather unconventional 
move since the momentum flux tensor does not originate 
from a microscopic collisional invariant. This move is 
motivated by the fact that the search for good entropies 
cannot keep up with the needs raised by various Lattice 
Boltzmann models: while it is rather easy to establish the 
set of constraints on the equilibrium, it is much less easy, 
even approximately, to find an entropy whose maximum, 
under fixed conserved fields, would also imply the rest of 
the constraints which do not come from the conservation 
laws. For this reason, one could try to force the solution, 
proceeding with any entropy of choice, but with more 
constraints. For instance, it is clear that minimizers of 
convex functions under constraints, 

J2fP[l,Vi,v i v i ]=p[l,u,(P/p)l + uu], (38) 

i 

are Galilean-invariant by construction since they encode 
the correct equilibrium form of the momentum-flux ten- 
sor right at the outset. (Any further requirements on 
the equilibrium can be added in the same manner.) The 
picture is clear: for D spatial dimensions, the above set 
of constraints generates iV c = 1 + D + D(D + l)/2 = 
(D + 1)(D + 2)/2 equations for the set of N c Lagrange 
multipliers A,B,C forming the quasi- invariant Qi = 
A + B ■ v,i + C : V{Vi, where the symbol ":" stands for 
the tensor scalar product. These equations arc gener- 
ally non-linear (if the entropy function is not quadratic 
in the populations) and consequently very hard to solve- 
if solvable at all- analytically in order to deliver closed 
expressions for the overconstrained local equilibrium as 
a function of the hydrodynamic fields. But even this is 
not the main drawback. More importantly, such over- 
constrained equilibria confine the domain of the phase 



space where the entropy production is positive to rather 
"thin" subsets (Karlin and Succi, 1998). This explains 
why the Lattice Boltzmann method offers less stability 
in such cases than non- isothermal hydrodynamics, where 
constructions of the equilibria have been focused mostly 
on satisfying the conservation constraints regardless of 
their origin. 



VI. DIRECTIONS FOR FUTURE RESEARCH 

To date, there arc still several directions where the 
Lattice Boltzmann theory calls for substantial progress 
and upgrades. Here, we shall briefly touch upon just two 
major streams of development: 

• Thermal flows 

• Non-ideal fluids 

In a nutshell, the problem is to guarantee a sufficient 
lattice symmetry to ensure the correct evolution of both 
kinetic and potential energy, without loosing numerical 
stability. This is fairly non-trivial, since it involves con- 
trol on higher-order local kinetic momenta, such as the 
heat flux, as well as on non-local information due to po- 
tential energy interactions. 

As a further prospective area of future development we 
mention 

• Quantum systems 

for which much less has been done to date. 



A. Thermal LBE's 

Thermal field theories are notoriously hard (Umezawa, 
1993), and Lattice Boltzmann theory is no exception. To 
date, isothermal, or better said, athermal LB systems 2 
are way better understood than their thermal counter- 
parts. The correct treatment of energy degrees of free- 
dom in a lattice still poses a number of basic questions. 

Energy dynamics in Lattice Boltzmann models is ac- 
counted for by enlarging the set of discrete speeds, 26 
being the number of kinetic moments to be matched 
by the set of discrete speeds. Early experiments (Mc- 
Namara, Garcia and Alder, 1995) showed that thermo- 
hydrodynamic LBE's exhibit a high degree of molecular 
individualism] errors in higher order moments seem to 
penetrate down into the low-order thermo-hydrodynamic 



2 Since LBE is based on a collection of mono-energetic beams 
S(v — Vi) in velocity space, the discrete LB distribution function 
is more appropriately classified as a zero-temperature system. 
We shall nonetheless stick to the more intuitive, if less rigorous, 
notation of isothermal schemes, to indicate that temperature is 
not a dynamic variable. 
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manifold. Similar type of difficulties have been recog- 
nized since long also in continuum kinetic theory (Lewis, 
1967). 

Despite significant progress, starting with the work of 
Alexander et al (Alexander, Chen and Sterling, 1993) 
thermal LBE's remain less robust than their isothermal 
counterparts. 

Thermal LBE's might be subject to a sort of numerical 
instability similar to the one affecting high-order finite 
difference schemes: due to the presence of high-energy 
particles, large-set of discrete speeds give rise to high- 
order dispersion relations, which in turn might develop 
spurious branches and unphysical solutions. 

This points even more strongly to the importance of 
a lattice H theorem. Again, a qualitative effect of lat- 
tice discreteness enters the scene. In the continuum, 
the local H theorem (positive entropy production due 
to collisions alone) does not explicitly involve transition 
probabilities from pre- to post-collisional states, and the 
detailed balance among different states are determined 
only by particle distributions with no other weighting 
factors. As a result, the H theorem is unaffected by 
particle advection since entropy is passively transported 
along collision-free trajectories. Consequently the local 
H theorem automatically implies the overall H theorem 
(or "global H theorem"). In the discrete case, however, 
the H theorem does depend explicitly on transition proba- 
bilities, which are functions of local collisional invariants 
(Molvig et al, 1988; Teixeira, 1992; Chen, Teixeira and 
Molvig, 1997; Chen, 1995; Chen, 1997; Chen and Teix- 
eira, 2000) . This is necessary to ensure that the resulting 
local equilibrium has a suitable form for producing cor- 
rect thermo-hydrodynamic equations satisfying Galilean 
invariance (Frisch, Hasslacher and Pomeau, 1986; Frisch 
et al., 1987). The condition for achieving the correct ther- 
mohydrodynamics requires the equilibrium form given by 
Chen (1995), namely an expansion up to 0(u 3 ) of a dis- 
crete local Maxwellian: 

C = p 9l (T)e^ 7 ^exp[-( Vl - u) 2 /2T] (39) 

where T is the local hydrodynamic temperature, 

p(DT + u 2 )/2 = J2^ 

i 

This form can be realized via explicit lattice gas par- 
ticle collisions involving transition probabilities of gi(T) 
(Chen, 1995; Chen, 1997). On the other hand, as per the 
discussion above, relaxation to a given local equilibrium 
needs not proceed through explicit collisions. It can be 
shown that ( |39"| ) admits a local H function of the form: 

h = Y,fMh/9i) (40) 

i 

The global H function may be simply defined as 

ff = J>(*). 

X 



Because of its explicit dependence on transition proba- 
bilities, the H value changes during the LBE advection 
phase, the result being that as particles hop in discrete 
steps from one lattice site to another, the values of the 
H function change ahead of the transition probabilities 
(which change during advection). This hinders transla- 
tional invariance and sets the stage for a spontaneous 
symmetry breaking of the H theorem with no counter- 
part in the continuum. Since the transition probabili- 
ties, gt, are functions of local hydrodynamic temperature, 
their value changes from place to place as particles move 
around the lattice. The violation of the global H theorem 
manifests itself whenever the temperature field acquires 
a spatial dependence. On the other hand, though the 
transition probabilities are not unity, they become con- 
stants in the isothermal LB models, which explains why 
the above problem disappears for isothermal LB models. 
This essential difference between thermal and isothermal 
lattice models lies at the heart of the significantly better 
stability of the latter. 

In summary, aside from low-Mach number expansions 
in the usual lattice BGK, isothermal models have been 
shown to possess an H theorem (i.e., global H theorem), 
while thermal models have not. Proving an H theorem 
and subsequently constructing an H theorem obeying the 
Lattice Boltzmann process for thermal models remains 
an outstanding problem in Lattice Boltzmann research. 



B. LBE for non-ideal fluids 

Lattice Boltzmann schemes for multiphase, multicom- 
ponent fluid flows are often heralded as a most promis- 
ing territory for LBE research. Indeed, there exist sev- 
eral extensions of the plain Lattice Boltzmann schemes 
for hydrodynamics which are able to include non-ideal 
gas effects (intermolecular interactions) (Shan and Chen, 
1993, 1994 Swift, Osborne and Yeomans, 1995, Luo 1998, 
2000). Most of these models can be cast in the form of 
a generalized LBE where the right-hand side is enriched 
with a (self-consistent) source term Ff. 

fi{x + v i ,t + l)-f i (x,t) = -u\f i -f^](x,t) + F i (41) 

Formally, Fi is the lattice version of the one-body self- 
consistent force (vector notation relaxed for simplicity): 

F(l)f(l) = f /i 2 (l,2)Fi 2 (l,2)d2 (42) 

where /12 is the two-body distribution function, F12 the 
two-body force, and 1, 2 stand for six-dimensional phase- 
space coordinates. Since the two-body distribution is 
likely to be computationally intractable, Lattice Boltz- 
mann research moved in the direction of a dynamic lat- 
tice density functional theory (Hansen and Mc Donald, 
1986), in which the effective force F(l) is represented by 
a semi-empirical functional of the macroscopic state of 
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the system, namely density, flow fields and their deriva- 
tives: 

F(l) = -V*[p,u,...] (43) 

These schemes have been validated for a series of test 
cases and complex flows applications, yielding fairly in- 
teresting results (as an illustration, see Figure 7). 

However, none of them has been proved to admit 
an underlying H theorem (Rothman, Keller and Gun- 
stensen, 1991; Shan and Chen, 1993; Shan and Chen, 
1994; Swift, Orlandini and Yeomans, 1995). This is in 
part due to the fact that equilibrium distribution func- 
tions depend on non-local properties. In other words, 
non-local information determines the transition proba- 
bilities. Consequently, similar to the problem encoun- 
tered in the thermal LB models, advection changes glob- 
ally defined quantities such as a H function. Another 
major theoretical issue is the nature of the hydrody- 
namic limit near thin-interfaces where the low-Knudsen 
assumption behind LBE is seriously challenged (a similar 
issue arises in the formulation of LBE-based renormaliza- 
tion group treatments of fluid turbulence (Chen, Succi, 
Orszag, 1999)). 

Finally, as indicated earlier in this paper, the proper- 
ties of the H theorem in the presence of a self-consistent 
force, are much less understood even in the continuum, 
Even though a global entropy production in the sense of 
decreasing H still exists in continuum Boltzmann, the 
issue of convexity and uniqueness of a minimum point 
becomes moot because the local condition dH / dt < no 
longer guarantees that collisional dynamics will attain 
a global minimum. Besides serving the purpose of for- 
mulating better Lattice Boltzmann models, we believe 
this can be regarded as a fundamental research topic of 
general interest in non-equilibrium statistical mechanics 
(Lebowitz, 1999). 

C. LBE for quantum systems 

To date, Lattice Boltzmann research is largely dom- 
inated by non-quantum physics applications. Nonethe- 
less, LBE can and has indeed been extended to (sim- 
ple) cases of quantum mechanical motion, typically in the 
form of the non-relativistic Schroedinger equation (Succi 
and Benzi, 1993, Succi 1996, Meyer 1997, Boghosian 
and Taylor 1997). The stepping stone to the quantum 
Lattice Boltzmann equation is the identification of the 
four-spinor wave function 'J, obeying the quantum rela- 
tivists Dirac equation with a complex discrete distribu- 
tion function. By writing the Dirac equation as a com- 
plex Lattice Boltzmann equation, it can be shown that 
the non-relativistic limit taking the Dirac equation into 
the Schroedinger equation is formally analogous to the 
adiabatic limit yielding hydrodynamics from the Lattice 
Boltzmann equation. 

The proof proceeds by writing the Dirac equation (in 
one dimension for simplicity) for a relativistic particle of 



rest mass m as two coupled LBE's for a pair of complex 
upward and downward propagating spinors u(z, t) and 
d(z, t) respectively: 

d t u + d z u = u) c d (44) 
d t d — d z d = —uu c u 

where uo c = mc 2 /h is the Compton frequency of a mate- 
rial particle. By applying the unitary transformation: 

cj> ± = -^(u±id) e*"' 

the above set of equations takes the following hydrody- 
namic form: 

d t (j> + + d z <t>- = (45) 
dt(j>~ — d z cf) + = 2iuo c u 

From these equations it is apparent that the "slow" (hy- 
drodynamic) symmetric mode is conserved, whereas the 
"fast" antisymmetric mode oscillates with a doubled fre- 
quency 2uj c . It is now readily checked that the adiabatic 
assumption: 

\dt<f>~\ « \2iw c <f>-\ 

delivers precisely the Schroedinger equation for a free- 
particle of mass to. Since the adiabatic approximation 
involves an imaginary frequency, unlike the classical case, 
the hydrodynamic branch (i.e., the Schroedinger equa- 
tion) remains time-reversible. This reversibility is some- 
how weaker than for the Dirac equation since it only ap- 
plies to real time, whereas the Dirac equation is reversible 
in both real and imaginary time. Interesting connections 
to a quantum if -theorem might arise in the framework 
of the ./V-body quantum Lattice Boltzmann equation, a 
totally unexplored issue to the best of our knowledge. 

The quantum LBE can be turned into a practical 
numerical scheme for the non-relativistic Schroedinger 
equation as well as for relativistic wave mechanics (Succi, 
2002). It might also be suitable to quantum computing 
paradigms. 

Saying that quantum LBE is a major issue of future 
Lattice Boltzmann research is probably an overstate- 
ment. Still, owing to the growing interest in quantum 
computing, the subject might be worth some attention 
in the years to come. 

VII. CONCLUSIONS 

These considerations complete the description of the 
Lattice Boltzmann method for isothermal Navier-Stokcs 
equations as a self-contained kinetic theory admitting a 
proper if theorem. The corresponding theory for fully 
thermo-hydrodynamic LBE's is still in its infancy. Not 
only is further understanding of this subject crucial to 
the formulation of more stable numerical algorithms, but 
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it might also stimulate new insights into fundamental 
physical questions involving temperature dynamics and 
non-local intermolecular interactions in discrete dynam- 
ical systems. 
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FIG. 1 The BBGKY hierarchy and its lattice analogue. Lack 
of microscopic detail becomes less and less relevant as one 
proceeds upwards along the hierarchy (from Succi, (2001)). 



FIG. 2 The hexagonal lattice of the Frisch-Hasslacher- 
Pomeau cellular automaton. Particles move along the six 
discrete links and meet at lattice nodes where they interact 
according to mass and momentum-conserving collision rules. 



FIG. 3 A typical collision in the Frisch-Hasslacher-Pomeau 
cellular automaton. Both post-collisional outcomes (right) 
are equally probable. The left hand hexagon shows two par- 
ticles initially at 9 o'clock and 3 o'clock moving in one time 
step to the central site. In the next time step, indicated by 
the hexagons on the right, their collision (conserving energy 
and momentum) results in the particles occupying the sites at 
11 o'clock and 5 o'clock or 1 o'clock and 7 o'clock. They could 
return to their initial locations at 9 o'clock and 3 o'clock but 
this event is ruled out since it would not produce any observ- 
able effect since the particles are indistinguishable. 



FIG. 4 Collisional relaxation procedure. Curves represent en- 
tropy levels, surrounding the local equilibrium / cq . The solid 
curve L is the entropy level with the value H(f) = H(f*), 
where / is the initial, and /* is the conjugate population. 
The vector A represents the collision integral, the sharp an- 
gle between A and the vector — VH reflects the entropy pro- 
duction inequality. The point M is the minimum entropy 
state on the segment [/,/*] The result of the collision up- 
date is represented by the point /(/3). The choice of /3 shown 
corresponds to the "overrelaxation" : H(f(/3)) > H(M) but 
H{f{P)) < H(f). The particular case of the BGK collision 
(not shown) would be represented by a vector Abgk, pointing 
from / towards / cq , in which case M = / cq . 
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FIG. 5 Evolution of a one-dimensional front in a shock tube. 
Density profile (dimensionlcss lattice units) is shown at t — 
500 for viscosity v — 3.3333 x 10~ 2 . The figure compares 
three Lattice Boltzmann algorithms on the lattice with 800 
nodes. Thin line: Exact solution at zero viscosity. Solid 
circles: Simulation. ELBM: The Lattice Boltzmann method 
based on the entropy function (from Ansumali and Karlin 
(2000)). LBGK: The Lattice BGK algorithm based on the 
polynomial ansatz of Qian, d'Humieres and Lallemand (1992). 
LBE: The Lattice Boltzmann model of Qian, d'Humieres and 
Lallemand (1991). The value of viscosity is taken close to the 
instability of the LBE. 




FIG. 6 Velocity profile in the one-dimensional shock tube 
benchmark. Simulation setup and notation are the same as 
in Figure 5 
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FIG. 7 A typical fluid interface exhibiting a single-mode 
Ray leigh- Taylor instability associated with a heavy fluid sit- 
ting on top of a light one. The density ratio of heavy and 
light fluids is 3. Re = (gL) 1/2 L/v = 1024, t/To = 3.5 
(To = Lf(gL) 1 ^ 2 ), g being the gravitational acceleration driv- 
ing the instability and L the box size. Periodic boundary con- 
ditions are used in horizontal directions and no slip boundary 
conditions are used in vertical directions (from Zhang, Chen 
and Doolen (1999), courtesy of Raoyang Zhang). 



